Introdução

O modelo de crescimento exponencial incorpora o princípio fundamental de que a vida provém de outra e assim examina o comportamento de uma população que pode crescer sem limitação alguma. Ao incluir a influência que o aumento na densidade populacional pode ter na limitação do crescimento, o modelo logístico incorpora outro princípio fundamental que é o da perda de excedente de indivíduos esperado por um crescimento ilimitado. Caso uma população crescesse de forma exponencial todo o universo conhecido poderia estar colonizado por um único tipo de vida. Outro componente essencial da modulação dos eventos demográficos é a interação com outras espécies.

Os modelos de competição de populações em denso dependência do arcabouço de reações de difusão (estocásticos) diferem dos modelos baseados em taxas constantes (determinísticos), pois os eventos demográficos são representados como eventos probabilísticos que podem ou não ocorrer em intervalos de tempo muito pequenos. Dessa forma, os modelos estocásticos possibilitam eventos de extinção de populações que poderiam permanecer indefinidamente no predito determinístico; ou flutuações nos tamanhos populacionais, propriedade ausente nos modelos determinísticos de denso dependência linear. Os modelos estocásticos podem ser descritos por uma equação mestra, um conjunto de equações diferenciais de primeira ordem que descrevem como a probabilidade de todos os possíveis estados do sistema varia no tempo contínuo. A partir da equação mestra é possível calcular a equação diferencial do primeiro momento da distribuição de probabilidades do sistema. Usando a aproximação de campo médio, onde é pressuposto variância zero, a equação diferencial do primeiro momento converge para a solução determinística. Alternativamente à descrição de todos os estados do sistema, é possível simular trajetórias individuais no sistema com o algoritmo de Gillespie (1977), onde podemos acompanhar a trajetória da abundância de cada espécie até a extinção.

Aqui vamos explorar a competição de 2 populações em crescimento denso dependente, N e M, em cenários de invasão da espécie M no sistema dominado por N que se encontra em sua capacidade de suporte. Usando modelos estocásticos e a aproximação de campo médio, explorei cenários de assimetria competitiva e de tamanho inicial da espécie invasora. Me pergunto a) Como a probabilidade de substituição da espécie residente pela invasora varia com o tamanho populacional inicial da invasora e da assimetria competitiva? e b) Como a presença da espécie residente reduz o tempo de extinção da espécie invasora?

Métodos

Reações

Vamos modelar a mudança de estado do sistema, descrito por (n,m), onde n e m são a abundância de N e M, respectivamente. O modelo considera possíveis mudanças discretas de 1 indivíduo no sistema em um intervalo de tempo muito pequeno dt. As oito reações consideradas estão disponíveis na tabela 1, elas constituem um modelo de denso dependência em que a competição interespecífica é modelada de forma similar à da competição intraespecífica: pela probabilidade de exclusão de um dos dois indivíduos por evento de colisão de dois indivíduos. O modelo pressupõe que não é possível distinguir entre indivíduos de uma mesma espécie, assim, o parâmetro \(\gamma\) descreve a probabilidade de qualquer um dos dois indivíduos ser excluído, por isso ele está dividido por 2.

Tabela 1. Possíveis reações no modelo estocástico de competição de duas espécies sob denso dependência
reações descrição P(reação)/dt unidade do parâmetro
n -> n + n reprod. assex. de N \(\beta n\) \(P(nasc. de N)/ind. de N dt\)
n -> 0 morte de N \(\delta n\) \(P(morte de N)/ind. de N dt\)
n + n -> n comp. intraesp. de N \(\gamma n (n-1) 2^{-1}\) \(P(exclusão de N)/colisão N dt\)
n + m -> m exclusão de N por M \(\text{C}_{n} n m\) \(P(exclusão de N)/colisão N M dt\)
m -> m + m reprod. assex. de M \(\beta m\) \(P(nasc. de M)/ind. de M dt\)
m -> 0 morte de M \(\delta m\) \(P(morte de M)/ind. de M dt\)
m + m -> m comp. intraesp. de M \(\gamma m (m-1) 2^{-1}\) \(P(exclusão de M)/colisão M dt\)
m + n -> n exclusão de M por N \(\text{C}_{m} m n\) \(P(exclusão de M)/colisão M N dt\)

The master equation and the mean-field approximation

A partir das reações é possível desenvolver uma equação mestra que descreve a taxa de mudança da probabilidade de um determinado estado [\(\frac{\partial P(n,m,t)}{\partial t}\)] pela diferença na taxa da probabilidade de entrada e saída deste estado (equação 1). Na primeira linha depois da igualdade, há a soma das taxa com que N pode perder um indivíduo, indo do estado (n+1,m) para o estado (n,m). Na segunda linha a taxa com que N pode ganhar um indivíduo, indo do estado (n-1,m) para o estado (n,m). Nas linhas 3 e 4 a mesma lógica para M. Na quinta linha há a taxa de saída do estado (n,m). Para detalhes veja tabela 1.

\[ \frac{\partial P(n,m,t)}{\partial t} = \\ P(n+1,m) (n+1) (\delta_n + n \gamma_n / 2 + m C_{n}) +\\ P(n-1,m) (n-1) \beta_n+\\ P(n,m+1) (m+1) (\delta_m + m \gamma_m / 2 + n C_{m}) +\\ P(n,m-1) (m-1) \beta_m+\\ - P(n,m)[n(\beta_n + \delta_n + (n-1) \gamma_n / 2 + m C_{n}) + m(\beta_m + \delta_m + (m-1) \gamma_m / 2 + n C_{m})] \]

Na aproximação de campo médio buscamos obter o primeiro momento da distribuição de probabilidade marginal de N, E[N] = , e de M, E[M] = . Para isso é necessário estabelecer algumas definições: \(<n> = \sum_{n=0}^{\infty} n P(n,m,t)\); \(<n^2> = \sum_{n=0}^{\infty} n^2 P(n,m,t)\); \(<m> = \sum_{m=0}^{\infty} m P(n,m,t)\); \(<m^2> = \sum_{m=0}^{\infty} m^2 P(n,m,t)\); \(<nm> = \sum_{n=0}^{\infty}\sum_{m=0}^{\infty} nm P(n,m,t)\). Para obter o primeiro momento marginal de N ou M basta multiplicar a equação mestra por n ou m e realizar o somatório (apêndice 1), e então é possível obter as equações 2a e 2b:

\[ \frac{d <n>}{d t} = <n> (\beta_n - \delta_n + \gamma_n2^{-1})- <n^2>\gamma_n2^{-1} - <nm>C_n \] \[ \frac{d <m>}{d t} = <m> (\beta_m - \delta_m + \gamma_m2^{-1})- <m^2>\gamma_m2^{-1} - <mn>C_m \]

A aproximação de campo médio pressupõe que \(<n>^2 = <n^2>\), \(<m>^2 = <m^2>\), \(<nm> = <n><m>\). Esse pressuposto implica que a variância da dinâmica populacional denso dependente é zero e que as populações são independentes. Outro pressuposto possível para simplificar a aproximação é de que \(\beta\) ou \(\delta\) >> \(\gamma\), assim \(\gamma\) pode ser omitido como fator do primeiro momento. Essa aproximação é útil pois permite a convergência das equações acima com o modelo determinístico de competição interespecífica de duas populações com crescimento denso dependente (eq. 3a e 3b). A investigação desse modelo determinístico mostra o quê poderíamos esperar no equilíbrio caso não houvesse nenhuma flutuação devido à estocástica dos eventos demográficos.

\[ \frac{d <n>}{d t} = <n> (\beta_n - \delta_n)- <n>^2\gamma_n2^{-1} - <n><m>C_n \] \[ \frac{d <m>}{d t} = <m> (\beta_m - \delta_m)- <m>^2\gamma_m2^{-1} - <m><n>C_m \]

The computational model

A ideia do algoritmo de Gillespie (1977) se baseia no fato de que a distribuição de tempos entre mudanças de estados pode ser descrita por uma distribuição exponencial, então computamos separadamente as mudanças de estado e o tempo entre mudanças, reduzindo o tempo de simulação. O tempo entre mudança de estados pode ser obtido pela divisão do log de uma amostra de distribuição uniforme padrão pela soma das taxas de saída do estado (janela de código 1). Existem 4 possíveis transições de estado: cada uma das duas espécies pode ganhar ou perder um indivíduo por evento demográfico no intervalo de tempo dt. A taxa de ganho de um indivíduo de N (\(n\_p\)) é \(n\_p=\beta n\). A taxa de perda de um indivíduo de N (\(n\_m\)) é \(n\_m = n[\delta +(n-1)\gamma2^{-1} + mC_n]\). A mesma lógica se aplica para as taxas de ganho e perda da espécie M (eq. 1 e janela de código 1). Após o sorteio do tempo do próximo evento demográfico, é feito um sorteio da mudança de estado com probabilidade dada pela contribuição relativa de cada taxa de saída (janela de código 1), então o sistema é atualizado. Nossa simulação termina quando alguma população atinge o zero.

janela de código 1

f_2sppSLV <- function(N0,beta_n,delta_n,gamma_n,C_n,
                      M0,beta_m,delta_m,gamma_m,C_m,
                      replicas,path_csv){
 df_ab <- data.frame(n=N0,m=M0,t=0)
 df_t <- data.frame(n=c(1,-1,0,0), # 4 possibles transitions:
                    m=c(0,0,1,-1)) # n + 1, n - 1, m + 1, m - 1
 f_loop <- function(){
   while(df_ab[nrow(df_ab),"n"]!=0 & df_ab[nrow(df_ab),"m"] !=0){
   v_rates <- c( 
     n_p = df_ab[nrow(df_ab),"n"] * beta_n,     # n+1
     n_m = df_ab[nrow(df_ab),"n"] * (delta_n +  # n-1
                                       (df_ab[nrow(df_ab),"n"]-1) * gamma_n/2 +
                                       C_n * df_ab[nrow(df_ab),"m"]),
     m_p = df_ab[nrow(df_ab),"m"] * beta_m,     # m+1
     m_m = df_ab[nrow(df_ab),"m"] *(delta_m +   # m-1
                                      (df_ab[nrow(df_ab),"m"]-1) * gamma_m/2 +  
                                      C_m * df_ab[nrow(df_ab),"n"]))
   v_Srates <- sum(v_rates)
   # update
   df_ab[nrow(df_ab)+1,] <- c(df_ab[nrow(df_ab),1:2] + df_t[sample(1:4,size = 1,prob=v_rates/v_Srates),],
                              df_ab[nrow(df_ab),"t"] - log(runif(1))/v_Srates)
   }
   return(df_ab)
}
df_dinamica <- plyr::rdply(.n = replicas,f_loop())
readr::write_csv(df_dinamica,file = path_csv)
}

Cenários Simulados: assimetria competitiva e tamanho inicial da invasora

Apenas o parâmetro que descreve a probabilidade de exclusão de uma espécie pela outra por colisão (parâmetros \(C_n\) e \(C_m\), tabela 1) variou entre simulações. Os demais parâmetros populacionais foram iguais: \(\beta =\) 1.916e-03, \(\delta =\) 1.889e-03, \(\gamma =\) 1.095e-06. Os parâmetros de competição interespecífica (\(C_n\) e \(C_m\)) variam entre \(\gamma\) e \(2\gamma\), tal que a assimetria competitiva, definido por \((C_n - C_m)/ \gamma\), variou entre -1 e 1 (figura 1). O tamanho inicial da espécie invasora variou entre 1% e 100% da capacidade de suporte K = 50 indivíduos (Figura 1). Foram simulados 11 cenários de tamanho inicial da invasora e 11 cenários de assimetria competitiva, totalizando 121 cenários, cada um com 100 réplicas.

Figura 1 Cenários simulados de competição interespecífica. Eixo y: assimetria competitiva - a diferença no parâmetro C da espécie residente pela invasora, dividido pelo parâmetro \(\gamma\). Eixo x: tamanhos iniciais da espécie invasora de 1% até 100%.

anaĺise dos dados

Para responder a primeira pergunta, como a probabilidade de substituição da residente varia em função do tamanho inicial da invasora e da assimetria competitiva, é possível usar a abordagem de campo médio e a de simulação de trajetórias individuais. No caso da abordagem de campo médio, exceto pelo cenário de equivalência das espécies (assimetria competitiva nula) que depende totalmente das condições iniciais e do regime de perturbações, dado tempo suficiente haverá certamente a substituição ou a não substituição da espécie residente pela invasora, em função da assimetria competitiva a favor da invasora ou da residente, respectivamente. Nessa abordagem basta plotar as soluções de equilíbrio das equações 3a e 3b em um plano de fase definido pelas abundância de N e M para avaliar qual o ponto de equilíbrio estável para definir se haverá ou não evento de substituição. No caso da simulação de trajetórias individuais, considerei a identidade da espécie que permanece após a primeira extinção. Se a espécie que permanece é a invasora então ocorreu um evento de substituição da residente pela invasora. Para descrever a probabilidade de substituição na simulação das trajetórias indivíduais usei um GLM binomial com o tamanho inicial da invasora e a assimetria competitiva como variáveis preditoras. Usei uma abordagem baseada em modelo médio (REF), em que a partir de um modelo cheio é ajustado todos os possíveis submodelos e então a predição de cada modelo é ponderada pelo peso de evidência do submodelo. O GLM binomial cheio considerou a interação do polinômio de segundo grau de cada preditora.

Para responder a segunda pergunta, como a presença da residente reduz o tempo de extinção da invasora, é necessário descrever o tempo de extinção da invasora nas simulações de trajetórias individuais em duas situações. A primeira situação é obtida pela análise da distribuição de tempos de extinção da invasora na presença da competidora, ou seja, os tempos totais até a espécie invasora se extinguir, seja ela a primeira ou última a se extinguir. Essa primeira situação pode ser descrita por um LMM em que os tempos de extinção de réplicas de uma mesma bateria de simulação (combinação única de tamanho inicial e assimetria competitiva) são agrupadas na estrutura aleatória (REF), e então usei as mesmas preditoras do GLM binomial para construir um modelo cheio. Na segunda situação considerei apenas o tempo parcial de extinção na ausência da residente, ou seja, considerei que o tempo é zerado no momento do evento de substituição. Para descrever a segunda situação é possível usar um LM em que o modelo cheio é composto pelo polinômio de segundo grau do tamanho da espécie invasora no momento do evento de substituição. Em ambas situações usei uma abordagem de modelo médio para construir o predito para os tempos de extinção. Então calculei a diferença no predito pelo modelo médio na segunda situação pela primeira situação.

Resultados

Como a probabilidade de substituição varia com o tamanho inicial da invasora e a assimetria competitiva?

PLANOS:

  • gráfico 1 plotar as isolinhas no plano de fase N X M em 3 situações extremas: C_nm = 1,0,-1, paras M0 = k*(0.01,0.25,0.5,0.75,1)

  • gráfico 2 predito para a Probabilidade de substituição em função de C

Figura 2 Gráfico Exploratório da Probabilidade da espécie invasora substituir a espécie residente em função da assimetria competitiva (eixo y) e tamanho populacional inicial da espécie invasora (eixo x).

Como a presença da residente reduz o tempo de extinção da invasora?

# Discussão

Uma expectativa é que a probabilidade da espécie invasora substituir a espécie residente deve aumentar quanto maior o tamanho populacional inicial da invasora e maior a assimetria competitiva em favor da invasora. Com a intenção de avaliar se é possível ter alguma intuição sobre o comportamento do sistema estocástico a partir da aproximação de campo médio, me pergunto como o predito determinístico diverge do observado em uma série de simulações individuais do processo estocástico. Pois, quanto menor o tamanho populacional maior a chance de uma sequência de eventos demográficos leve à extinção da invasora, mesmo quando dominante. E quanto maior a assimetria competitiva em favor da invasora, maior a probabilidade da espécie residente reduzir em abundância por possível evento demográfico.

Apêndice

Modelos estatísticos

Probabilidade de substituição da residente pela invasora

Figura A1 Proporção observada de substituições

O modelo cheio mais plausível é aquele que considera a interação entre os polinômios (Tabela 1). Porêm este modelo não apresenta um excelente ajuste aos dados (Figura SI 1 e sumário do módelo 1). O predito e observado diferem em valores extremos: quando a proporção de substituições é baixa, o predito pode subestimar a probabilidade de substituição (Figura SI 2); quando a proporção é alta, o predito pode superestimar a probabilidade. Para as análises utilizei os pacotes DHARMa (REF), MuMIn (REF) e AICcmodavg (REF).

tabela A1 seleção de GLMs Binomiais para descrever a probabilidade de substituição

Figura A2 Resíduos Quantílicos GLM Binomial cheio. Resíduos quantílícos obtidos com a função simulateResiduals do pacote DHARMa, com 250 simulações.

tabela A2 Sumário do modelo 1: GLM Binomial Cheio

## 
## Call:
## glm(formula = cbind(c_invasion, 100 - c_invasion) ~ (C_nm.z + 
##     I(C_nm.z^2)) * (M0.z + I(M0.z^2)), family = "binomial", data = df_resultados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2851  -0.7129  -0.1248   0.5908   2.7772  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -0.69649    0.04619 -15.079  < 2e-16 ***
## C_nm.z                 0.21493    0.03138   6.850 7.39e-12 ***
## I(C_nm.z^2)           -0.03052    0.03547  -0.860  0.38966    
## M0.z                   0.86425    0.03928  22.000  < 2e-16 ***
## I(M0.z^2)             -0.26621    0.04088  -6.512 7.40e-11 ***
## C_nm.z:M0.z            0.07671    0.02821   2.719  0.00655 ** 
## C_nm.z:I(M0.z^2)       0.01055    0.02870   0.368  0.71316    
## I(C_nm.z^2):M0.z       0.04878    0.03156   1.545  0.12225    
## I(C_nm.z^2):I(M0.z^2) -0.04809    0.03226  -1.491  0.13607    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1842.58  on 120  degrees of freedom
## Residual deviance:  139.94  on 112  degrees of freedom
## AIC: 711.01
## 
## Number of Fisher Scoring iterations: 4

Modelo Médio

Figura A3 predito pelo modelo médio e observado.

Tempo de Extinção

Tempo para extinção na ausência de competidora

Figura A4 Conjunto de dados completo do Tempo estimado de extinção na ausência da competidora

Ajuste modelo cheio LM TE 1sp

Tabela A6 Tabela seleção do modelo referência para descrever o log(Tempo de Extinção).

##           dAICc  df weight
## log(M0)^2    0.0 4  1     
## log(M0)     25.1 3  <0.001
## M0^2       197.0 4  <0.001
## M0         548.6 3  <0.001
## 1         1675.2 2  <0.001

Figura A5 Resíduos Quantílicos LM referência para descrever o tempo de extinção na ausência de residente. Resíduos quantílícos obtidos com a função simulateResiduals do pacote DHARMa, com 250 simulações.

tabela A6 Sumário do modelo 2: LM referência polinômio de segundo grau

## 
## Call:
## lm(formula = log_TE ~ log_M0 + I(log_M0^2), data = df_1sp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0319 -0.7538 -0.0595  0.7187  4.3239 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.33084    0.14619  43.306  < 2e-16 ***
## log_M0       1.43874    0.10589  13.588  < 2e-16 ***
## I(log_M0^2) -0.09455    0.01814  -5.213 1.97e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.058 on 3351 degrees of freedom
## Multiple R-squared:  0.3939, Adjusted R-squared:  0.3935 
## F-statistic:  1089 on 2 and 3351 DF,  p-value: < 2.2e-16

Modelo Médio

Ajuste modelo cheio GLM Gamma(log) TE 1sp

Tabela A7 GLM Gamma(log): Tabela seleção do modelo referência para descrever o log(Tempo de Extinção).

##             dAICc  df weight
## N:log(M0)^2    0.0 4  1     
## N:M0^2        67.5 4  <0.001
## log(M0)^2    159.8 4  <0.001
## log(M0)      186.8 3  <0.001
## M0^2         257.9 4  <0.001
## M0           519.8 3  <0.001
## 1           1348.4 2  <0.001

Figura A5 Resíduos Quantílicos LM referência para descrever o tempo de extinção na ausência de residente. Resíduos quantílícos obtidos com a função simulateResiduals do pacote DHARMa, com 250 simulações.

tabela A6 Sumário do modelo 2: LM referência polinômio de segundo grau

## 
## Call:
## glm(formula = log_TE ~ log_M0 + I(log_M0^2), family = Gamma(log), 
##     data = df_1sp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.49937  -0.08442  -0.00735   0.07504   0.49797  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.872745   0.018798  99.623  < 2e-16 ***
## log_M0       0.182179   0.014996  12.148  < 2e-16 ***
## I(log_M0^2) -0.015193   0.002816  -5.396  7.4e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.01418678)
## 
##     Null deviance: 59.818  on 2746  degrees of freedom
## Residual deviance: 38.801  on 2744  degrees of freedom
## AIC: 8524.6
## 
## Number of Fisher Scoring iterations: 4

Figura A6 gráficos diagnostico do lm: log(TE) ~ log(M0)^2

Modelo Médio

Figura A7 predito e observado pelo modelo médio: lm - log(TE) ~ log(M0)^2

Predição final Modelo Médio LM: log(TE) ~ log(M0)^2

Figura A8 Predito pelo modelo médio (lm: log(TE) ~ log(M0)) e o observado na escala log(TE) e na escala padrão TE. Em vermelho a média e em azul o intervalo de confiança computado pressupondo normalidade assintótica (REF Modavgpred)

Tempo de extinção na presença da residente

Figura A9 Gráficos exploratórios do Tempo de Extinção total: quando a substituição foi um sucesso ou fracasso

Tabela A7 Seleção do modelo cheio para descrever o tempo de extinção total.

##                    dAICc df weight
## C_nm^2 + log(M0)^2   0.0 7  1     
## C_nm^2 * log(M0)^2  26.6 11 <0.001
## log(M0)^2           59.9 5  <0.001
## C_nm^2 + M0^2      200.4 7  <0.001
## M0^2               203.4 5  <0.001
## C_nm^2 * M0^2      228.5 11 <0.001
## 1                  486.6 3  <0.001
## C_nm^2             495.3 5  <0.001

Figura A10 Resíduos quantílicos lmer TEI 2spp cheio: C_nm^2 + log(M0)^2

Tabela A8 sumário do lmer TEI 2spp cheio mais plausível: C_nm^2 + log(M0)^2

## Linear mixed model fit by REML ['lmerMod']
## Formula: log_TE ~ (C_nm.z + I(C_nm.z^2)) + (log_M0.z + I(log_M0.z^2)) +  
##     (1 | id)
##    Data: df_2sp
## 
## REML criterion at convergence: 35042.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2272 -0.6977 -0.0942  0.5977  4.9156 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.000    0.000   
##  Residual             1.057    1.028   
## Number of obs: 12100, groups:  id, 121
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    9.432693   0.017075 552.429
## C_nm.z         0.087216   0.009347   9.331
## I(C_nm.z^2)   -0.049695   0.010584  -4.695
## log_M0.z       0.797119   0.013809  57.725
## I(log_M0.z^2) -0.076225   0.009602  -7.938
## 
## Correlation of Fixed Effects:
##             (Intr) C_nm.z I(C_.^ lg_M0.
## C_nm.z       0.000                     
## I(C_nm.z^2) -0.620  0.000              
## log_M0.z    -0.414  0.000  0.000       
## I(lg_M0.^2) -0.562  0.000  0.000  0.736
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')